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We extended the Density Matrix Renormalization Group method to study the real time dynamics 
of interacting one dimensional spinless Fermi systems by applying the full time evolution operator to 
an initial state. As an example we describe the propagation of a density excitation in an interacting 
clean system and the transport through an interacting nano structure. 



I. INTRODUCTION 

The density matrix renormalization group method 
(DMRG)^ is a powerful technique to study the prop- 
erties of one dimensional interacting quantum systems. 
The advantage of the DMRG is that it can treat quantum 
lattice systems in the presence of site-dependent inter- 
action, hopping parameter, and on-site potentials^ with 
high accuracy, including subtle lattice effects like multi- 
ple umklapp processes^ 

Originally the method was set up to describe the equi- 
librium properties of the ground state and a few ex- 
cited states. It was then extended to calculate fre- 
quency dependent spectral functions by the use of the 
[H - E Q - lu + irj\ ~ 1 operator 

A second approach to study transport properties 
within the framework of DMRG is to relate equilibrium 
properties of the ground state to transport properties. 
Molina et al& and Meden and SchollwockS calculated 
the conductance through an interacting nano structure 
attached to leads by relating the conductance of the sys- 
tem to the ground state curvature, based on an idea by 
Sushkovia 

Cazalilla and Marstonii used the basis states of the 
last DMRG sweeps to integrate the Schrodingcr equa- 
tion in real time. However, as shown by Luo, Xiang, 
and Wangjia their approach was inconsistent with the 
DMRG scheme, since they used a basis to integrate the 
Schrodingcr equation which was only adapted to the 
t = state, neglecting the relevant states for t > during 
the DMRG sweeps. 



II. CALCULATION OF THE TIME EVOLUTION 

Instead of integrating the time dependent Schrodinger 
equation numerically, we make use of the formal solution 
and apply the full time evolution operator 

U(t 2 ,h) = e -^-^ (1) 

to calculate the time dependence of an initial state |£(0)) 



|e(t)>=U(t > 0)|f(0)), (2) 
where IH is the Hamiltonian of the system of interest. 



While the calculation of e _i:K * is not feasible for large 
matrix dimensions, one can calculate the action of a ma- 
trix exponential on a vector similar to the diagonaliza- 
tion of sparse matrices, where one cannot diagonalizc 
the full matrix, but one can search for selected eigen- 
values and eigenvectors. We apply a Krylov subspace 
approximation^ to calculate the action of the matrix ex- 
ponential in equation Q on a state £). We would like to 
encourage the reader interested in implementing a ma- 
trix exponential to study the excellent review by Moler 
and Van Loanii and discourage the use of a Taylor ex- 
pansion. In our implementation we make use of a Pade 
approximation from Expokit— to calculate the dense 
matrix exponential in the Krylov subspace. In order to 
calculate |£(t)) up to a final time T, we discretize the 
time interval into N time steps {tQ,t\,t2, ■■ ■ , i/v} with 
<o = 0, tpf = T and tj < ij+i, typically tj — tj-i = 0.5 
and N = 2T. It turns out that using a time slice tj — tj-i 
of the order of one is sufficient to ensure a fast conver- 
gence of 



\^t j ))=e-^- t ^\at j -i)) 



(3) 



It is crucial that one does not need to store the wave func- 
tion at all time steps. Instead one can add imme- 
diately to the density matrix p and calculate the matrix 
elements of interest at each time step separately. There- 
fore, only the two wave functions and are 
needed during the calculation of the time evolution^ 

One even does not need to calculate the ground state 
of the unperturbed system, it is sufficient to keep the ma- 
trix elements of the 3~C which may be useful for systems 
with highly degenerate ground states, e.g., an XXZ fer- 
romagnet, where the convergence of the diagonalization 
is very slow. However, in this work we calculate the m 
lowest lying eigenstates \^ m ) of JC, E m \^ m ) = < K\ i & m ), 
and use the time evolution operator only on the subspace 
orthogonal to these eigenstates, 



(4) 



m=0 
m— 1 



m) = £ e-(^-*)*|* ra ><* w |«0)> (5) 

m=0 

+ e -«K-Eo)t(l-p)\Z(0)) l 

and calculate the dynamics in the subspace of the eigen- 
states {^0: ■ ■ • j ^m-i} exactly. In addition we intro- 
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duced a phase choice e lEot to make the ground state time 
independent. 

III. WAVE PACKET DYNAMICS 

In order to prepare an initial state we apply a small 
perturbation <5IH to the Hamiltonian !H of the system 
of interest and calculate |£(0)) as the ground state of 
'K + S'K. 

As a first example we study the time evolution of a 
density pulse in a model of interacting spinless fermions: 

M M 

m = -* Y c x°x-i + 4-i c z + V Y n * n x-i ( 6 ) 

X — 1 X— 1 

where t is the hopping parameter and V the nearest 
neighbour interaction parameter. In this work we mea- 
sure all energies with respect to t = 1. To create a wave 
packet we add a Gaussian potential 

M , 

E («-«!) 
e ^ (7) 

where \i is the strength, a the width, and the posi- 
tion of the perturbation. In Fig. we have plotted the 
time evolution of an initial wave packet at X\ = 6 in a 
30 site system at half filling and periodic boundary con- 
dition using 300 states per DMRG block. Due to time 
reversal symmetry the initial state consist of a left and 
a right moving wave packet. During the time evolution 
the initial peak splits into two peaks which are moving 
with the group velocities ±v g . For V — the DMRG 
results coincides with the result from exact diagonaliza- 
tion. As a first result this method gives direct access 
to the group velocity v g of a density excitation without 
relying on finite size analysis or arguments from confor- 
mal field theory>i£ In table {1} we compare the extracted 
group velocities v g for the M = 30 site system with the 
Fermi velocity vf known from Bcthc ansatz results for an 
infinitesimal excitation in the infinite system size limit, 
Vf = 7rsin(2ry)/(7r — 2r]) with the usual parametrization 
V = — 2 cos(2f7)pil As expected from vp the wave pack- 
ets travel faster the stronger repulsive interaction are, 
while they are slowed down by attractive interaction. For 
fi = 0.002 there is a good agreement from v g with vf, 
while the results for /i = 0.02 already include dispersion 
effects. In addition, the broadening of the wave packets 
reveals information on the dispersion relation. A detailed 
study is beyond the scope of this work and subject for 
future studies. 

It is not obvious that one can target for a few low ly- 
ing states | 'I'm) of "K, the ground state £(0) of IK + <55f 
and N ~ 100 time steps of simultaneously in each 
DMRG step. However, since the DMRG truncation is 
the only approximation in our method, we can systemat- 
ically increase the number of states n n c u t kept per block 
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TABLE I: Comparison of v g extracted from DMRG simu- 
lations for a M = 30 site system and a potential strength 
fi = 0.02 , /j, = 0.002 and Bethe ansatz results for vf in the 
infinite system and an infinitesimal small excitation. 

to control errors due to the Hilbert space truncation. In 
Fig. J2J we plot an initial wave packet £(0) and the wave 
packet £(T) at T = 8 for a 50 site system, a potential 
strength of // = 0.02, an interaction strength V — 1.0 
and periodic boundary conditions. While for the initial 
state 200 states per block seem to be sufficient to describe 
the wave packet, far more states are needed to obtain the 
dynamics of the wave packet correctly. The slow conver- 
gence at the boundaries of the system is related to an 
implementation detail. We do not keep all density oper- 
ators to evaluate n x in the final iteration step. Instead 
we calculate n x during the last 1.5 finite lattice sweeps, 
when we have the operators available. The advantage of 
this procedure is that one does not need to keep all oper- 
ators at the price that the operators close to boundaries 
have less accuracy, since they are evaluated in a highly 
asymmetric block configuration. 

Remarkably, the overlap (^o|£(0)) between the ground 
state of "K and the ground state of the system of 3i + 
83< shown in Fig. @ is 99.99%. Therefore, it is the 
0.01% contribution which gives the initial excitation and 
governs the time evolution. This high overlap was the 
motivation to introduce the projection defined im Eq. 10}. 

IV. TRANSPORT THROUGH A QUANTUM 
DOT 

In order to study transport through an interacting 
nano structure, we prepare a system consisting of an 
interacting region coupled to nonintcracting leads, see 
Fig. ©, 

M 712-1 

ft = -* Yl ( c x c x-i + c l-i c x) + U n xn x -i 

X — 1 X—7li-\-l 

+ lU Y n x n x -r, (8) 

x—ni ,n-2 

where t is hopping parameter, U is the interaction on the 
nano structure and 7 defines a smoothening of the onset 
of interaction at the nano structure. In this work we have 
set 7 = 0.5, compare Molina et al&. In the following we 
denote Mg = ni — n\ the number of sites in the nano 
structure, M the number of site of the total system and 
Ml = M — Ms the number of lead sites. 

In Fig. 10} we show the time evolution of a wave packet 
initially placed in the left lead of a system with Ms = 7, 
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FIG. 1: Time evolution of a wave packet for a system of M = 30 sites, periodic boundary conditions, and V = 0.0 (plus), 
V = 1.0 (crosses) and V = 2.0 (stars). The snapshots are taken at T — 0, 2, 5, and 8. ncut = 300 states per block were kept 
within the DMRG procedure. 
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FIG. 2: (a) Initial wave packet for a 50 site system with pe- 
riodic boundary conditions and V = 1, for different numbers 
of states kept per block: nc u t = 100 (plus), 200 (crosses), 300 
(stars), 400 (open squares), 500 (filled squares) and 750 (open 
circles), (b) Same system at time T — 8. 



ing sites to smoothen out the Friedel oscillations (for all 
V) and the charge density wave on the dot for V = 5.0. 
At the beginning of the time evolution, the wave packet 
is not overlapping with the interacting nano structure, 
hence the packets travel synchronously for all interaction 
strength. Once they reach the nano structure, the wave 
packets move with the group velocity of the interacting 
system. After the wave packets have left the interact- 
ing region they continue to move at the velocity of the 
noninteracting system. The wave packets which traveled 
through the interacting region are now travelling in front 
of those with smaller interaction. 

For U <= 2.0 the nano structure is transparent, al- 
though there seems to be a reflection of a negative pulse 
as predicted by Safi and Schulai^. For very strong in- 
teraction, U > 2.0, there is an instability to a charge 
density ordering^ii and the nano structure has a finite 
reflection. We would like to remark that these simula- 
tions clearly demonstrate that the Luttinger description 
of the infinite system makes already sense for a system 
consisting of a few lattice sites only. One should should 
keep in mind that for such small systems the effective 
parameters, like vf, have not reached the infinite system 
limit. However, the scaling already leaves its fingerprint. 



M = 50, interaction V = 0.0, 1.0, 2.0 and 5.0 and hard 
wall boundary conditions, which lead to perfect reflec- 
tion at the chain ends. To rule out truncation errors we 
use ncut = 1000 states per block, compare discussion of 
Fig. J2J). We have averaged the density over neighbour- 



V. SUMMARY 

In summary we have shown an accurate method to 
calculate real time dynamics within the framework of 
DMRG. By applying the matrix exponential on our ini- 
tial state state we are able to perform the complete time 
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FIG. 3: Nano structure attached to leads. 

integration of the time dependent Schrodinger equation 
in each single DMRG step. In this setup the only ap- 
proximation is given by the truncation procedure of the 
DMRG, which be can systematically checked by increas- 
ing the number of states kept during the DMRG sweeps. 

VI. NOTE 

While preparing this work we became aware of 
related workiSi 2 ^ on using real time dynamics within 



the DMRG. Both apply a Suzuki- Trotter decomposition 
of the time evolution operator, based on an work by 
by Vidali^i In addition, their work relies on the state 
prediction^ to calculate the time evolution of a state, 
which represents an additional approximation. In 
contrast to their work we calculate the initial state and 
apply the full time evolution operator in each iteration 
step, without introducing additional approximations 
beyond the truncation scheme of the DMRG. 
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FIG. 4: Transport through an interacting region of Ms — 7 sites and Ml = 43 lead sites, hard wall boundary conditions and 
n-Cut = 1000. The snapshots are taken at T = 0, 3, 7, 12. The density is averaged over 2 neighbouring site to smoothen 2fc_F 
oscillations and is plotted vs. the site location x for V = (plus), V = 1 (crosses), V — 2 stars, and V = 5 (squares). For 
V — the line is calculated by an exact diagonalization. 



